回帰係数の多重比較

単回帰 y=a+bx における，H₁:a=0 と H₂:b=0 の多重比較

data(thuesen, package = "ISwR")
d <- thuesen
str(d)

## 'data.frame':    24 obs. of  2 variables:
##  $ blood.glucose : num  15.3 10.8 8.1 19.5 7.2 5.3 9.3 11.1 7.5 12.2 ...
##  $ short.velocity: num  1.76 1.34 1.27 1.47 1.27 1.49 1.31 1.09 1.18 1.22 ...

res <- lm(short.velocity ~ blood.glucose, d)
ctab <- as.data.frame(summary(res)$coef)
round(ctab, 4)

##               Estimate Std. Error t value Pr(>|t|)
## (Intercept)      1.098     0.1175   9.345   0.0000
## blood.glucose    0.022     0.0105   2.101   0.0479

b <- coef(res)
V <- vcov(res)
m <- nrow(d) - length(b)
C <- diag(2)
s <- sqrt(diag(C %*% V %*% t(C)))
t <- C %*% b/s
R <- cov2cor(C %*% V %*% t(C))

# q : adjusted p value
library(mvtnorm)
q <- sapply(abs(t), function(x) 1 - pmvt(-c(x, x), c(x, x), corr = R, df = m))
ctab$adjust.p <- q
round(ctab, 4)

##               Estimate Std. Error t value Pr(>|t|) adjust.p
## (Intercept)      1.098     0.1175   9.345   0.0000   0.0000
## blood.glucose    0.022     0.0105   2.101   0.0479   0.0637

# u : critical value for siginificant level 0.05
u <- qmvt(p = 0.95, tail = "both.tails", df = m, corr = R, delta = c(0, 0))$quantile
# CI for 95% family-wise confidence level
ctab$lwr <- C %*% b - u * s
ctab$upr <- C %*% b + u * s
round(ctab, 4)

##               Estimate Std. Error t value Pr(>|t|) adjust.p     lwr    upr
## (Intercept)      1.098     0.1175   9.345   0.0000   0.0000  0.8367 1.3590
## blood.glucose    0.022     0.0105   2.101   0.0479   0.0637 -0.0013 0.0452